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Modeling Dynamically Coupled Fluid-Duct Systems 
with Finite Line Elements 


J.B. Saxon 

Rockwell International, Space Systems Division 
Huntsville, Alabama 35806 


Abstract 

Structural analysis of piping systems, especially dynamic analysis, typically 
considers the duct structure and the contained fluid column separately. Coupling 
of these two systems, however, forms a new dynamic system with characteristics 
not necessarily described by the superposition of the two component system's 
characteristics. Methods for modeling the two coupled components simultaneously 
using finite line elements are presented. Techniques for general duct intersections, 
area or direction changes, long radius bends, hydraulic losses, and hydraulic 
impedances are discussed. An example problem and results involving time 
transients are presented. Additionally, a program to enhance post-processing of 
line element models is discussed. 

Introduction 

Structural analyses of piping systems are usually accomplished with finite 
line element models to which estimated fluid loads are applied. Under transient 
conditions, however, dynamic coupling often occurs between the pipe and the 
contained fluid such that the two can no longer be considered separately. In the 
analyses of Space Shuttle Main Engine propellant feedlines and ducts, the need 
for modeling structural dynamic response of coupled fluid-duct systems has led to 
new analysis methods. These methods allow simultaneous analysis of the 
structure and contained fluid column by representing both as overlaying strings of 
line elements. This requires the fluid to be treated as one dimensional; an 
assumption considered adequate from the structural analyst's perspective. The 
same basic approach has been previously explored (Reference 1), but the 
specifics detailed herein were developed independently and contain some unique 
features which represent an expanded application of the underlying principles. 
The bulk of the presented methods deals with maintaining the correct force transfer 
between the duct and fluid elements. The theory includes treatment of forces 
transferred at general intersections or direction/area changes, forces transferred at 
long radius bends, forces due to losses, and terminal hydraulic impedances. Post- 
processing animations are also discussed. 

Dynamic Approach 

In the presence of both steady-state and transient dynamics, it is desirable to 
consider the two separately and superimpose their results for a complete answer. 
In the proposed method, the steady-state dynamics are analyzed with a quasi-static 
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approach for which the FEA representation of the fluid column assumes static 
equilibrium, but is understood to represent some steady-state dynamic condition. 
Forces transferred by the quasi-static model accurately mirror those of the true 
steady-state dynamic condition except that head losses are usually ignored. If a 
loss is considered significant, then the steady-state analysis should include forces 
applied by the analyst to the fluid column and/or structure as a correction. Losses 

will be addressed later in greater detail. 

The underlying assumption of the transient analysis is that transient dynamic 
loads produce the same response whether applied to a quasi-static system or to a 
genuinely steady-state dynamic system. Under transient loads, velocities and 
loads experienced by the FEA representation of the fluid are understood to add 
vectorially to the corresponding steady-state values. Forces due to losses are 
generally a function of velocity squared, and therefore should have a transient 
component as well. Discussion of a method for estimating transient losses will also 
be deferred, but will be built on the assumption that transient velocities are small 

compared to steady-state velocity. ... 

By analyzing steady-state and transients separately, a quantifiable error is 
introduced in the rate at which pressure perturbances are propagated up and 
downstream. In reality, pressure wave propagation should be a superposition of 
acoustic and steady-state velocities. Error due to the absence of a steady-state 
velocity is dismissed, however, since acoustic velocities dwarf flow velocities 
within the scope considered here. 

Modeling a Fluid Column and Coupled Straight Pipe 

Analysis of a duct or pipe-line structure with finite line elements is 
commonplace. The fluid coupling methods presented herein are intended as an 
addition to the standard structural model of a duct. Therefore, definition of the 
structural portion of the coupled system requires no particular deviation from 
standard practice, nor will it be discussed in any detail. 

Treatment of the contained fluid column is one dimensional and therefore 
representable by structural rod elements (structural members that carry only 
tension/compression). Properties of the rod elements are based on properties of 
the fluid column they represent. Mass density of the rod should equal density of the 
fluid, and cross-sectional area of the rod should equal the cross-sectional area of 
fluid'column. Young's Modulus of the rod corresponds to Effective Bulk Modulus, 
Be, of the fluid, which is calculated from Bulk Modulus, B, corrected for radial 
stiffness of the pipe as shown below. 

Be = (BEt)/(2BR + Et) 0> 

where E = Young's Modulus of pipe 
t = pipe wall thickness 
R = nominal radius of pipe 

For duct cross-sections more complex than cylindrical pipe, such as bellows. Be 
may be calculated by methods such as demonstrated in Reference 2. With these 
fluid/structural analogies in place, the axial stress in the rod is taken to equal the 
total pressure, static plus dynamic, in the fluid column. 

The validity of this approach can be confirmed by a simple model of a fluid 
column constrained in all but the axial direction along its length and completely 
constrained at one end. FEA modal analysis should match hand calculated open- 
closed organ pipe frequencies. 
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For the sake of conveniently coupling fluid to pipe, it is beneficial to adopt 
certain conventions in defining the two coincident strings of line elements which 
represent the duct structure and the fluid column. First, it is convenient that nodes 
on one element string have its nodes coincident with nodes of the other string, such 
that pairs of corresponding fluid/structural nodes exist. Each node pair should be 
coupled such that they move independently in the direction corresponding to flow, 
but are otherwise rigidly joined. This is illustrated schematically in Figure 1, and 
may be accomplished with either Multi-Point-Constraint (MPC) equations or zero 
length springs. As a convenience to defining the coupling, nodal degrees of 
freedom, especially those of the fluid node, should be aligned with one axis in the 
direction of flow. As a matter of convention, it will be assumed herein that the nodal 
Z axis is so oriented. 



Figure 1. FEA of Straight Duct and Coupled Fluid Column 

Coupling at Direction/Area Changes 


At any point where the duct contains a bend, an intersection, or a change in 
flow area, considerations beyond those described thus far must be applied in order 
to correctly account for the transfer of forces between the fluid column and the duct 
structure at that location. For now, the discussion will be limited by the assumption 
that the duct feature in question can be represented as a point along the duct path. 
If the size of the feature is considerable compared with the overall model, namely, 
the case of a very long radius elbow, an alternate approach described later may be 
more appropriate. 

Figure 2 illustrates a pipe intersection involving all the features under 
consideration. (Planar geometry is not necessary, but is used here for clarity.) 
Mass continuity between the three fluid columns and for^e transfer between the 
fluid and the duct can be imposed with a single MPC equation involving the Z 
translations of each fluid node and the translational degrees of freedom of the 
intersection’s structural node. The equation is derived by treating the intersection's 
structural node as if it were another incoming fluid branch oriented and sized such 
that all the Pressure-Area forces acting on the intersection added vectorially to 
zero. 

A generalized approach, as applied to Figure 2, proceeds as follows. Each 
of the three branches entering the intersection should be modeled as described 
previously, and each branch should have its own fluid node at the point of 
intersection. For convenience, direct the Z axes of the intersection's three fluid 
nodes into the intersection. Let the vectors Aj equal the flow area of branch i 
times the unit Z vector of branch i's fluid node at the intersection. Each Aj should 
be expressed in the reference frame defining the nodal displacements of the 
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intersection's structural node. Now define a vector As according to 

3 

A S + I(Aj) = 0 

i=1 


Having defined A s , the MPC equation enforcing mass continuity and force balance 


is 
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(Asx’Xs) + (Asy’Ys) + (Asz*Zs) + X(Aj*Zj) = 0 

i=i 


(3) 


STRUCTURAL NODE 
AT INTERSECTION 



Figure 2. Duct Intersection and Area Change Example 


Of course, equations 2 and 3 are applicable to any number of incoming 
branches. If applied to one branch, a capped end is defined. If applied to two 
branches, an elbow is defined. 

For some analysis software, it may be necessary to enforce these constraints 
with large stiffness values rather than an MPC equation. In effect, the stiffness 
equations form the intersection element defined below for n branches intersecting 
at a point. Here, Ki is an arbitrary multiplier just large enough to make fluid at the 
point of intersection relatively incompressible. 
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Long Radius Elbows 

in some cases, the radius of a bend is not negligible compared to the total 
length of duct being modeled. Figure 3 shows such a bend, encompassing an 

angle 0, as modeled with n fluid/structural node pairs spaced evenly between the 
end points of the bend. As with straight pipe, the fluid and structural nodes should 
move independently in the direction tangent to the direction of flow, but be rigidly 
coupled perpendicular to flow as illustrated. This arrangement accurately models 
the forces transferred, however, the stresses in the rod elements will exceed the 
total pressure in the fluid. This can be corrected by adjusting the properties of the 
rods so as to maintain correct axial stiffness and achieve correct stresses. The 

basic properties of the rods in the bend, e.g. A, E, and p, should be adjusted to A’, 
E', and p' according to 

A'=A/sin(0/(n+1)) ( 5 ) 

E'=E«sin(0/(n+1)) ( 6 ) 

p'=p*sin(0/(n+1)) 


Figure 3. FEA of Curved Duct and Coupled Fluid Column 
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Losses 


In addition to force transfer between fluid and duct at direction/area changes, 
head losses may play a significant role in how forces are distributed throughout the 
system. As previously stated, steady-state losses must be accounted for by the 
analyst . This is accomplished with equal and opposite forces applied to the fluid 
column in the upstream direction, and to the structure in the downstream direction. 
The magnitude of the loss, where steady-state velocity is a given quantity, must be 
determined as shown below, where C is obtained from reference data. 

FLoss=A*C»p»(Vss)^ ^ 


Figure 4 illustrates an elbow under steady-state conditions with an assumed loss. 
The mass and force balance equations would make pressure in the fluid column 
the same both up- and downstream of the elbow. The additionally applied forces 
for loss allow the fluid column downstream to see a reduced pressure while 
maintaining mass continuity. 



In the presence of transients, the total velocity, Vss+Vt, should produce a 
Force of loss equal to the sum of the steady-state loss, Fss. and a transient 
component of loss, Ft, as shown below. 

Floss - F ss + F t =A-C-p-(V ss+ v,)2 < 9 > 

Assuming that the loss variable C does not appreciably change between velocities 
Vss and Vss+Vt, substituting equation 8 into 9 and solving for Ft gives 

Ft = A-Op-(2*V SS ‘Vt+Vt 2 ) (10) 
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If a further assumption is made that Vt 2 is insignificant compared with 2*Vss*Vt, the 
transient force of loss can be estimated as 

F t = (2A*Op*V ss )*Vt t 11 ) 

which is a linear damping term. As long as C is constant, the above estimation is 
90% accurate for Vt < (.2*Vss)- The damper implied above should be affixed at the 
point of loss between the structural node and the downstream fluid node as seen in 
Figure 4. The damper increases the force of loss for a transient velocity that 
increases flow velocity magnitude and decreases the force of loss for a transient 
velocity that decreases flow velocity magnitude. Unfortunately, the above treatment 
does not account for the fact that much of the loss represented by published values 
of C accounts for irreversible losses several diameters downstream of the elbow, 
thus reducing the static pressure in the fluid without causing an equal and opposite 
load on the duct. None-the-less, the above method does provide an estimate of 
force distribution due to transient losses. 


Terminal Hydraulic Impedances 


The ducts addressed to date have all terminated at the inlet to a rocket 
engine or an accumulator device which provided a quantifiable impedance to flow. 
The Hydraulic impedance is expressed as compliance, resistance, and inertance, 
for which the mechanical equivalents are stiffness, resistance, and mass, 
respectively. For the steady state analysis, the terminal end of the fluid column can 
be fixed, thereby precluding the need to account for the impedance. In the trans- 
ient analysis, however, the fluid column must terminate at a spring/mass/damper 
system, as shown in Figure 5, with K, M, and D defined in terms of compliance, 
mass, and resistance as follows: 


K = r»A/Compliance 

(12) 

M = r»A 2 *lnertance 

(13) 

D= I> A 2 * Resistance 

(14) 


Where r is weight of density of fluid 


Figure 5. 
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Mass-Spring-Damper System to Model Terminal Hydraullic Impedance 


Note that these impedance terms are often frequency dependent. This requires the 
analyst to choose coefficients carefully, possibly analyzing frequency domain 
transients piece-wise over the frequency range of interest. 


Example Problem 


As a demonstration of the modeling concept, a transient analysis is 
presented based on work performed for NASA's Marshall Space Flight Center. 
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The Center is home to Technology Test Bed, a facility for test firing Space Shuttle 
Main Engines. The problem addresses the possibility of incorporating an 
oscillating piston in a side branch of the main Liquid Oxygen feedline in order to 
agitate pressure at the engine inlet (See Figure 6.) This test procedure, known as 
pogo pulsing, simulates the effect of space craft vibrations on the propellent 
systems, which, in turn, has been observed to cause and couple with engine thrust 
oscillations. The immediate problem was to predict pogo pulsing's effect on the 
facility feedline. 
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Figures 6. Schematic of TTB Liquid Oxygen Feedline 


The coupled fluid-duct model of the feedline incorporates many of the 
modeling features described herein. At both elbows and the pulser tee, stiffness 
matrices representing general area/direction changes were applied. Losses at 
valves, flowmeters, etc. were modeled with linear dampers, the damping values 
based partially on pressure data recorded during operation. The terminal 
hydraulic compliance of the engine was modeled with a spring, but, for the sake of 
this example, its frequency dependence and any accompanying hydraulic 
resistance were ignored. 

Quasi-static analysis of steady-state operation was relatively straightforward. 
The applied loads, including those to account for losses, are shown in Figure 7. 
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i TANK PRESSURE, 6319 Ibf 
¥ 



Figure 7. Loss Loads Applied to Feedline Model for Quasi-static Analysis of 

Steady-state Conditions 

The results of three Modal analyses, shown in Table 1, demonstrate the 
inadequacy of analyzing the fluid and structural system separately . The first set of 
modes considers the fluid acoustics only, as if the duct were perfectly rigid (except 
for radial expansion accounted for by equation 1). The second set of modes 
considers duct flexibility, but treats the fluid only as added mass. The last set of 
modes, which accounts for the fluid-duct coupling, is clearly more than the simple 
superposition of the first two cases. 

Modes from the coupled system were used in two transient analyses 
simulating the start-up of pogo pulser operation at 15 Hz and 35 Hz. Both analyses 
were driven by enforcing sinusoidal displacement of the piston node. Figures 8 
and 9 show response of the system presented as nodal displacement of a point in 
the fluid column near the engine inlet, and displacement of a node on the duct near 
the downstream elbow. The transient response to 15 Hz Pogo pulsing quickly 
settles to a steady-state containing higher frequency components - as might be 
expected. As might also be expected, response to 35 Hz pogo pulsing exhibits 
resonance with the coupled system mode at 34.9 Hz, amplifying the input signal. 

Obviously, any analysis of Pogo pulsing based on the separate models of 
the fluid and duct would have been much more crude. The fact that 35 Hz was a 
critical system frequency would have remained unknown since it shows up in 
neither of the first two columns of Table 1. Even at non-critical frequencies, the 
model can be used to quantify structural stresses as a function of the oscillating 
pressure seen at the engine inlet and/or the stroke of the piston. Hopefully, this 
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method will prove applicable and useful in the structural analysis of other coupled 
fluid-duct systems to address problems involving waterhammer, hydraulic control, 
accumulator systems, etc. Already, the method has been employed in the analysis 
of other NASA related systems such as the duct and accumulator, components of 
the space shuttle main engine, seen in Figures 12 and 13. 


Table 1. Comparison of Modal Results 


Fluid 

Acoustic Modes, 
Rigid Duct (Hz) 

Duct Structural 
Modes, 

-Frozen" Fluid (Hz) 

Coupled 
Fluid-Duct 
System Modes (Hz) 

1.8 


1.8 


6.9 

6.1 

29.6 


14.3 


32.4 

22.1 

59.1 


26.6 


66.1 

32.2 


77.4 1 

34.9 

88.5 


62.5 


92.6 

71.8 


97.6 

84.8 

118.0 


89.8 


124.2 

93.6 


133.6 

98.3 


142.2 

118.6 

147.6 


130.9 

177.0 


143.8 


181.1 

145.7 


192.9 

162.7 

205.2 


179.4 

I 


191.5 


Response to 15 Hz Pogo Pulsing 



Time (sec) 


fluid response structural response 


Figure 8. Transient Feedline Response to Pogo Pulsing Start-up at 15 Hz 
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Response to 35 Hz Pogo Pulsing 



Time (sec) 


Fluid Response Structural Response 


Figure 9. Transient Feedline Response to Pogo Pulsing Start up at 35 Hz 

Post-Processing Animations 

As is common to all line element models, deformed geometry animations of 
models such as described herein are bland. The information sought, namely 
pressure distribution in the fluid and its correlation with structural response, is 
difficult to convey visually with line element plots. In order to enhance 
understanding of the solution, a Fortran program called DUCT6D was written 
which, in conjunction with PATRAN 2.5, simultaneously presents the model, its 
structural deformations, and fluid pressures in a color animation. DUCT6D 
converts the line elements into rings of shells elements so that the circular cross 
section of the duct can be visually perceived. The color of the shells changes with 
the pressure in the fluid column. Deformation of the shell model mirrors the 
extrapolated displacements and rotations of the line model. DUCT6D will produce 
both modal and transient animations. DUCT6D is available from COSMIC (See 
NASA Tech Briefs . August 1 993, p. 33) in a format compatible with output from the 
EAL analysis code, however, minor editing of input formats should adapt it to other 
specific analysis packages. Figure 10 shows the undeformed Lox feedline from the 
example problem as depicted with DUCT6D. Figure 11 shows the 35 Hz mode 
shape from the example problem as depicted with the benefit of DUCT6D, where 
the darker color indicates higher pressure. Figures 12 and 13 further demonstrate 
DUCT6D’s usefulness in visualizing line element models. Figure 12 is a plot of a 
coupled fluid-duct model of a duct/accumulator system found on the space shuttle 
main engine. Figure 13 is the same model as enhanced with DUCT6D. 
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Figures 10. Feedline Model as 
Depicted with DUCT6D 



Figures 1 1. Mode 7 of Feedline Model, 
34.9 Hz, as Depicted with DUCT6D 



Figures 12. Line Element SSME 
Lox Duct 


Figures 13. SSME Lox Duct 

After DUCT6D 
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